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Abstract 


A  study  of  the  neutrally-stratified  flow  within  and  over  an  array  of  three-dimensional  (3-D) 
buildings  (cubes)  was  undertaken  using  simple  Reynolds-averaged  Navier-Stokes  (RANS) 
flow  models.  These  models  consist  of  a  general  solution  of  the  ensemble- aver  aged,  steady- 
state,  three-dimensional  Navier-Stokes  equations,  where  the  k-e  turbulence  model  (A:  is 
turbulence  kinetic  energy  and  e  is  viscous  dissipation  rate)  has  been  used  to  close  the 
system  of  equations.  Two  turbulence  closure  models  were  tested;  namely,  the  standard  and 
Kato-Launder  k~e  models.  The  latter  model  is  a  modified  k-e  model  designed  specifically  to 
overcome  the  stagnation  point  anomaly  in  flows  past  a  bluff  body  where  the  standard  k-e 
model  overpredicts  the  production  of  turbulence  kinetic  energy  near  the  stagnation  point. 
Results  of  a  detailed  comparison  between  a  wind  tunnel  experiment  and  the  RANS  flow 
model  predictions  are  presented.  More  specifically,  vertical  profiles  of  the  predicted  mean 
streamwise  velocity,  mean  vertical  velocity,  and  turbulence  kinetic  energy  at  a  number  of 
streamwise  locations  that  extend  from  the  impingement  zone  upstream  of  the  array,  through 
the  array  interior,  to  the  exit  region  downstream  of  the  array  are  presented  and  compared 
to  those  measured  in  the  wind  tunnel  experiment.  Generally,  the  numerical  predictions 
show  good  agreement  for  the  mean  flow  velocities.  The  turbulence  kinetic  energy  was 
underestimated  by  the  two  different  closure  models.  After  validation,  the  results  of  the 
high-resolution  RANS  flow  model  predictions  were  used  to  diagnose  the  dispersive  stress, 
within  and  above  the  building  array.  The  importance  of  dispersive  stresses,  which  arise  from 
point-to-point  variations  in  the  mean  flow  field,  relative  to  the  spatially-averaged  Reynolds 
stresses  are  assessed  for  the  building  array. 
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Resume 


Une  etude  des  eeoulements  d'air  a  stratification  neutre  a  l’interieur  et  par-dessus  une  ma- 
trice  de  batiments  (cubes)  tridimensionnelle  (3-D)  a  ete  entreprise  a  Faide  de  modeles 
d’ecoulements  simples  Navier-Stokes  de  Reynolds  moyens  (RANS).  Ces  modeles  consistent 
en  une  solution  generate  de  la  moyenne  de  Fensemble  des  equations  Navier-Stokes  tridimen- 
sionnelles,  en  etat  d’equilibre,  ou  le  schema  de  turbulence  k-e  ( k  est  1  energie  cinetique  de 
la  turbulence  et  e  est  le  taux  de  dissipation  visqueuse)  a  ete  utilise  pour  conclure  le  systeme 
des  equations.  Deux  schemas  de  fermeture  de  turbulences  ont  ete  testes:  le  schema  stan¬ 
dard  et  le  schema  k-e  Kato-Launder.  Ce  dernier  schema  est  un  schema  k-e  modifie,  congu 
specifiquement  pour  surmonter  Fanomalie  du  point  d’ arret  des  eeoulements  a  contour ner 
les  corps  non  profiles  ou  le  schema  standard  k-e  prevoit  des  valeurs  superieures  aux  chiffres 
observes  de  la  production  de  Fenergie  cinetique  de  la  turbulence  pres  du  point  d  arret. 
Des  resultats  d’une  comparaison  detaillee  entre  un  essai  de  la  soufflerie  aerodynamique 
et  les  predictions  du  modele  d'ecoulement  RANS  sont  presentes.  Plus  exactement,  les 
profils  verticaux  de  la  velocite  longitudinale  moyenne,  de  la  velocite  verticale  moyenne 
et  de  Fenergie  cinetique  de  la  turbulence,  prevus  a  un  certain  nombre  dtendroits  longi- 
tudinaux  qui  s?etendent  a  partir  de  la  zone  d  impaction  du  secteur  en  am  ont,  a  travers 
la  matriee  interieure,  jusqu’a  la  region  de  la  sortie  en  aval  de  la  matrice,  sont  presentees 
et  eomparees  a  celles  mesurees  dans  Fessai  de  la  soufflerie  aerodynamique.  Normalement, 
les  previsions  numeriques  correspondent  bien  avec  les  velocites  moyennes  des  eeoulements. 
L 5 energie  cinetique  de  la  turbulence  avait  ete  sous-estimee  par  deux  schemas  differents  de 
fermeture.  Apres  validation,  les  resultats  des  previsions  du  modele  d’ecoulement  RANS  de 
haute  resolution  ont  ete  utilises  pour  diagnostiquer  la  contrainte  dispersive,  a  l’interieur  et 
au-dessus  de  la  matriee  de  batiments,  L  importance  de  ces  eontraintes  dispersives  qui  se 
produisent  a  partir  de  variations  de  points  a  points  dans  le  champ  d  ecoulement  moyen, 
relatives  a  la  moyenne  spatiale  des  nombres  Reynolds  sont  evaluees  pour  la  matriee  de 
batiments. 
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Executive  summary 


Introduction:  It  is  anticipated  that  Canadian  Forces  (CF)  in  the  foreseeable  future  will 
have  to  fight  in  or  protect  urban  areas,  whether  in  battle,  peace-making,  peacekeeping, 
or  counter-terrorist  operations.  The  increased  awareness  and  importance  accorded  by  the 
public  worldwide  and  their  governments  to  maintain  appropriate  defences  against  chemical 
and  biological  warfare  (CBW)  agents  in  an  urban  (built-up)  environment,  the  prediction 
of  casualties  and  human  performance  degradation  resulting  from  such  releases,  and  the 
development  of  operational  procedures  and  regulations  to  control,  mitigate,  and  monitor  the 
fate  of  CBW  agents  in  urban  areas  with  high  population  densities,  will  require  mathematical 
modeling  of  urban  wind  flows  and  dispersion.  In  this  regard,  it  should  be  noted  that  the 
prediction  of  flows  in  the  urban  environment  is  in  principle  pre-requisite  to  or  co-requisite 
with  the  prediction  of  contaminant  (e.g..  CBW  agent)  dispersion  within  a  cityscape. 

Results:  A  study  of  the  neutrally-stratified  flow  in  and  over  an  array  of  three-dimensional 
(3-D)  buildings  (cubes)  was  undertaken  using  Reynolds-averaged  Navier-Stokes  (RANS) 
flow  models.  These  models  consist  of  a  general  solution  of  the  ensemble-averaged,  steady- 
state,  three-dimensional  Navier-Stokes  equations,  where  the  k-e  turbulence  model  ( k  is 
turbulence  kinetic  energy  and  e  is  viscous  dissipation  rate)  has  been  used  to  close  the 
system  of  equations.  Two  turbulence  closure  models  were  tested;  namely,  the  standard  and 
Kato-Launder  k-e  models.  The  latter  model  is  a  modified  k-e  model  designed  specifically  to 
overcome  the  stagnation  point  anomaly  in  flows  past  a  bluff  body  where  the  standard  k-e 
model  overpredicts  the  production  of  turbulence  kinetic  energy  near  the  stagnation  point. 
Results  of  a  detailed  comparison  between  a  wind  tunnel  experiment  and  the  RANS  flow 
model  predictions  are  presented.  More  specifically,  vertical  profiles  of  the  predicted  mean 
streamwise  velocity,  mean  vertical  velocity,  and  turbulence  kinetic  energy  at  a  number  of 
streamwise  locations  that  extend  from  the  impingement  zone  upstream  of  the  array,  through 
the  array  interior,  to  the  exit  region  downstream  of  the  array  are  presented  and  compared 
to  those  measured  in  the  wind  tunnel  experiment.  Generally,  the  numerical  predictions 
show  good  agreement  for  the  mean  flow  velocities.  The  turbulence  kinetic  energy  was 
underestimated  by  the  two  different  closure  models.  After  validation,  the  results  of  the 
high-resolution  RANS  flow  model  predictions  were  used  to  diagnose  the  dispersive  stress, 
within  and  above  the  building  array.  The  importance  of  dispersive  stresses,  which  arise  from 
point-to-point  variations  in  the  mean  flow  field,  relative  to  the  spatially-averaged  Reynolds 
stresses  are  assessed  for  the  building  array. 

Significance  and  Future  Plans:  The  k-e  turbulence  closure  model  with  an  isotropic, 
linear  eddy  viscosity  is  the  simplest  complete  turbulence  model  (viz.,  no  advance  knowl¬ 
edge  of  any  property  of  the  turbulence  is  required  for  the  simulation  other  than  the  initial 
and/ or  boundary  conditions  for  the  problem)  that  is  available  currently,  and  its  moderately 
good  predictive  performance  of  the  developing  flow  through  an  array  of  three-dimensional 
buildings  without  the  need  to  adjust  any  closure  constants  is  encouraging.  This  model 
may  be  useful  as  a  general-purpose  simulator  of  urban  flows  since  it  is  simple  enough  to  be 
tractable  numerically  and,  hence,  not  require  excessive  computing  time.  However,  before 
this  model  (or,  variants  of  it)  can- be  used  for  this  purpose,  it  needs  to  be  validated  against 
more  experimental  data  for  flow  in  various  obstacle  arrays.  From  this  perspective,  the 
availability  of  extensive  and  accurate  experimental  data  for  flow  through  various  obstacle 
arrays  is  urgently  required.  Finally,  the  accurate  prediction  of  urban  dispersion  requires 
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knowledge  of  the  mean  wind  and  turbulence  through  the  urban  area.  The  utility  of  the 
urban  flow  simulation  models  investigated  here  for  provision  of  the  disturbed  wind  field 
statistics  required  for  a  physically-based  dispersion  model  needs  further  investigation. 


Yee  E  and  Lien,  F.S.  (2004).  Numerical  Simulation  of  the  Disturbed  Flow  Through  a 
Three-Dimensional  Building  Array.  (DRDC  Suffield  TR  2004-108).  Defence  RtD  Canada 
-  Suffield. 


IV 


DRDC  Suffield  TR  2004-108 


Sommaire 


Introduction:  On  prevoit  que  les  Forces  canadiennes  (FC),  auront,  dans  un  avenir  assez 
rapproche,  a  combattre  a  Pinterieur  de  zones  urbaines  ou  a  proteger  ces  dernieres,  du- 
rant  des  operations  de  combat,  de  maintien  de  la  paix  ou  antiterroristes.  Une  prise  de 
conscience  accrue  ainsi  que  Fimportance  que  le  monde  entier  et  ses  gouvernements  accor- 
dent  au  maintien  de  moyens  de  defense  appropries  contre  les  agents  de  guerre  chimiques 
et  biologiques  (CB)  dans  les  milieux  urbains  (construits),  a  la  prevision  des  blesses  et  la 
degradation  de  la  performance  humaine  resultant  de  telles  emissions  ainsi  qu’a  la  mise  au 
point  de  procedures  operationnelles  et  reglements  de  controle  pour  attenuer  et  surveiller  le 
sort  des  agents  CB  dans  les  zones  urbaines  comprenant  des  hautes  densites  de  population, 
exigeront  des  modelisations  mathematiques  des  ecoulements  eoliens  urbains  et  de  leur  dis¬ 
persion.  A  cet  egard,  il  faut  noter  que  la  prevision  des  ecoulements  dans  un  milieu  urbain 
est  en  principe  pre-requise  ou  co-requise  avec  la  prevision  de  la  dispersion  du  contaminant 
(par  ex.:  agent  CB),  dans  un  paysage  urbain. 


Resultats:  Une  etude  des  ecoulements  d’air  a  stratification  neutre  a  Pinterieur  et  par-dessus 
une  matrice  de  batiments  (cubes)  tridimensionnelle  (3-D)  a  ete  entreprise  a  Paide  de  modeles 
d’ecoulements  simples  Navier-Stokes  de  Reynolds  moyens  (RANS).  Ces  modeles  consistent 
en  une  solution  generate  de  la  moyenne  de  Pensemble  des  equations  Navier-Stokes  tridimen- 
sionnelles,  en  etat  d’equilibre,  ou  le  schema  de  turbulence  k-e  ( k  est  Penergie  cinetique  de  la 
turbulence  et  e  est  le  taux  de  la  dissipation  visqueuse)  a  ete  utilise  pour  conclure  le  systeme 
des  equations.  Deux  schemas  de  fermeture  de  turbulences  ont  ete  testes  :  le  schema  stan¬ 
dard  et  le  schema  k-e  Kato-Launder.  Ce  dernier  schema  est  un  schema  k-e  modifie,  congu 
specifiquement  pour  surmonter  Panomalie  du  point  d’arret  des  ecoulements  a  contourner 
les  corps  non  profiles  ou  le  schema  standard  k-e  prevoit  des  valeurs  superieures  aux  chiffres 
observes  de  la  production  de  Penergie  cinetique  de  la  turbulence  pres  du  point  d’arret. 
Des  resultats  d’une  comparaison  detaillee  entre  un  essai  de  la  soufRerie  aerodynamique 
et  les  predictions  du  modele  d’ecoulement  RANS  sont  presentes.  Plus  exactement,  les 
profils  verticaux  de  la  velocite  longitudinale  moyenne,  de  la  velocite  verticale  moyenne 
et  de  Penergie  cinetique  de  la  turbulence  prevus  a  un  certain  nombre  d’endroits  longitu- 
dinaux  qui  s'etendent  a  partir  de  la  zone  d ’impaction  du  secteur  en  amont,  a  travers  la 
matrice  interieure,  jusqu’a  la  region  de  la  sortie  en  aval  de  la  matrice,  sont  presentees  et 
comparees  a  celles  mesurees  dans  Pessai  de  la  soufflerie  aerodynamique.  Normalement. 
les  previsions  numeriques  correspondent  bien  avec  les  velocites  moyennes  des  ecoulements. 
L’energie  cinetique  de  la  turbulence  aveit  ete  sous-estimee  par  deux  schemas  differents  de 
fermeture.  Apres  validation,  les  resultats  des  previsions  du  modele  d’ecoulement  RANS  de 
haute  resolution  ont  ete  utilises  pour  diagnostiquer  la  contrainte  dispersive,  a  Pinterieur  et 
au-dessus  de  la  matrice  de  batiments.  L’importance  de  ces  contraintes  dispersives  qui  se 
produisent  a  partir  de  variations  de  points  a  points  dans  le  champ  d’ecoulement  moyen, 
relatives  a  la  moyenne  spatiale  des  nombres  Reynolds  sont  evaluees  pour  la  matrice  de 
batiments. 

La  Portee  des  Resultats  et  les  Plans  Futurs:  Le  schema  de  fermeture  de  la  turbulence 
k-e  avec  une  viscosite  turbulente  lineaire  isotrope  est  le  modele  de  turbulence  complet  le 
plus  simple  (ziv.,  n’exigeant  aucune  connaissance  prealable  des  proprietes  de  la  turbulence 
pour  la  simulation  autre  que  celle  des  conditions  aux  limites  et/ou  initiales  du  probleme)  qui 
est  actuellement  disponible;  de  plus,  sa  performance  de  prevision  moderement  bonne  en  ce 
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qui  concerne  le  developpement  de  1 ’ecoulement  a  t ravers  la  matrice  des  batiments  tridimen- 
sionnels,  sans  avoir  besoin  d’ajuster  aucune  constante  de  fermeture,  est  encourageante,  Ce 
modele  pourrait  etre  utile  comme  simulateur  polyvalent  des  ecoulements  urbains  puisqu  il 
est  assez  simple  pour  etre  tractable  numeriquement  et  n' exige  pas  ainsi  de  temps  excessif  de 
ealcuL  Cependant,  avant  que  ce  modle  (ou  une  de  ses  variantes)  puisse  etre  utilise  dans  ce 
but,  il  doit  etre  valide  par  des  donnees  experiment  ales  supplementaires  pour  1  ecoulement 
dans  des  matrices  variees  d'obstacles.  Dans  cette  perspective,  il  est  urgent  de  trouver  des 
donnees  experiment  ales  etendues  et  exactes  pour  1  ecoulement  a  t  ravers  des  matrices  variees 
d’obstacles.  Enfin,  la  prevision  exacte  de  la  dispersion  en  milieu  urbain  exige  une  connais- 
sance  du  vent  moyen  et  de  la  turbulence  a  travers  ce  milieu  urbain.  L’utilite  des  modeles  de 
simulation  de  Pecoulement  en  milieu  urbain  examines  ici  pour  alimenter  les  statistiques  des 
champs  de  vent  perturbe  requis  pour  le  modele  de  dispersion  d’apres  les  eriteres  physiques, 
exige  d’etre  examinee  plus  longuement. 


Yee,  E.  and  Lien,  F.S.  (2004).  Numerical  Simulation  of  the  Disturbed  Flow  Through  a 
Three-Dimensional  Building  Array.  (DRDC  SufReld  TR  2004-108).  Defence  R&D  Canada 
-  Suffield. 
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Introduction 


The  surface  of  the  Earth  is  covered  with  roughness  elements  that  form  a  patchwork  mosaic 
of  varying  surface  roughness  such  as  crops,  forests,  and  urban  areas.  This  wide  range  of 
complex  surfaces  disturbs  the  turbulent  flow  over  the  surface  and  influences  the  processes 
that  govern  the  exchange  of  momentum,  heat,  and  mass  between  the  ‘complex'  surface  and 
the  atmosphere.  In  recent  years,  demographic  evolution,  cultural  and  technological  devel¬ 
opments,  and  economic  activity  have  accelerated  the  phenomenon  of  city  formation.  Huge 
urban  centres  that  consist  of  agglomerations  of  buildings,  often  with  irregular  geometry 
and  spacing,  present  some  of  the  most  ‘complex’  surfaces  for  which  an  understanding  of  the 
surface-atmosphere  interaction  and  the  nature  of  the  flow  between  the  obstacles  themselves 
are  not  currently  available. 


An  urban  canopy  consists  of  a  large  collection  of  buildings  and  other  obstacles  (e.g.,  cars 
lining  a  street,  treed  areas  in  city  green  spaces,  etc.)  that  are  aggregated  into  complex 
structures.  When  this  rough  surface  interacts  with  the  atmospheric  flow  within  find  above 
it,  the  disturbed  flow  field  can  become  extremely  complex.  Curved  mean  streamlines,  large 
velocity  gradients,  sharp  velocity  discontinuities,  flow  separations  and  reattachments,  cavity 
regions,  recirculation  zones,  and  strongly  inhomogeneous  turbulence  can  result  depending 
on  the  obstacle  geometry  and  spacing.  The  variety  of  mean  and  turbulent  flow  regimes  that 
result  from  the  flow  within  and  over  an  arbitrary  arrangement  of  non-uniform  obstacles 
is  too  large  to  be  catalogued.  It  would  be  naive  to  expect  a  universal  pattern  to  exist 
in  these  complex  disturbed  flows  and,  as  a  consequence,  a  comprehensive  physically-based 
theory  for  these  flows  does  not  exist  currently.  Nevertheless,  the  prediction  of  turbulent 
flows  through  an  urban  canopy  is  required  for  a  number  of  important  applications  (e.g., 
the  estimation  of  wind  loads  on  urban  structures  and  the  modelling  of  the  dispersion  of 
pollutants  in  or  near  an  urban  area).  Given  the  complexity  of  the  flow  and  turbulence 
structure  over  an  urban  surface,  where  flow  and  turbulence  lengthscales  and  timescales  can 
change  dramatically  from  one  obstacle  configuration  to  the  next,  it  seems  reasonable  to 
suggest  that  given  increases  in  both  computing  speed  and  memory  of  affordable  computers 
that  numerical  simulations  of  these  flows  might  provide  valuable  insights  that  could  make 
the  study  of  these  intractable  flows  easier  and  more  effective. 

Computational  fluid  dynamics  (CFD)  has  been  applied  to  the  problem  of  predicting  tur¬ 
bulent  flows  around  an  isolated  building  or  obstacle.  To  this  purpose,  two  methods  have 
been  used  to  predict  the  complex  flow  around  a  single  obstacle;  namely,  large-eddv  simu¬ 
lation  (LES)  in  which  the  large  eddies  are  computed  and  the  smallest,  subgrid-scale  eddies 
are  modelled  and  Reynolds-averaged  Navier-Stokes  (RANS)  equations  in  which  the  flow  is 
modelled  through  suitable  averaged  quantities  for  both  the  mean  and  turbulent  motion.  Pa¬ 
terson  and  Apelt  ([1],  [2])  used  a  RANS  model  (with  fc-e  turbulence  closure)  to  calculate  the 
mean  velocity  and  pressure  in  the  flow  around  a  single  prismatic  building  on  a  flat  surface. 
Zhang  et  al.  [3]  studied  the  effects  of  wind  shear  and  turbulence  in  the  upstream  flow  on 
the  development  of  the  flow  field  around  an  isolated  cuboid  (building)  using  the  standard 
k-e  model.  Yamada  [4]  applied  RANS  modelling  (using  a  higher-order  turbulence  closure) 
to  simulate  the  flow  around  a  single  square-based  building  on  complex  terrain.  Ferzinger 
and  Paric  [5]  describe  the  application  of  LES,  using  the  mixed  Smagorinsky  scale-similarity 
model  as  the  subgrid-scale  model,  for  an  investigation  of  the  unsteady  flow  over  a  cube 
mounted  on  one  wall  of  a  channel. 
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More  recently,  CFD  has  been  applied  to  simulate  the  flow  within  and  over  clusters  of 
buildings.  Certainly,  a  true  building  cluster,  made  up  of  a  group  of  buildings  of  roughly 
comparable  size,  is  expected  to  provide  the  most  complicated  flow  problem  for  CFD  because 
mutual  flow  interferences  between  buildings  need  to  be  taken  into  account  [6].  Much  of  the 
early  work  in  this  area  focussed  on  the  application  of  RANS  modelling  to  the  investigation  of 
the  three-dimensional  characteristics  of  flow  regimes  within  an  urban  canyon  (e.g.,  Dawson 
et  al.,  [7]:  Hunter  et  al,,  [8];  Baik  and  Kim,  [9]).  More  recently,  Liu  and  Barth  [10]  and 
Walton  and  Cheng  [11]  have  applied  LES  to  study  flow  and  pollutant  dispersion  in  an 
idealized  street  canyon.  Hanna  et  al.  [12]  used  LES  with  the  Smagorinsky  subgrid  closure 
model  to  simulate  the  mean  flow  and  turbulence  within  some  simple  obstacle  arrays.  Cheng 
et  al.  [13]  applied  both  RANS  (k-e  model)  and  LES  with  a  dynamic  Smagorinsky  subgrid 
model  to  predict  the  fully-developed  flow  over  a  matrix  of  cubes  immersed  in  a  plane  channel 
flow.  Lien  et  al.  [14]  applied  a  k-e  turbulence  closure  model  in  conjunction  with  linear  and 
nonlinear  eddy- viscosity  formulations  to  predict  the  developing  flow  over  an  array  of  two- 
dimensional  (2-D)  buildings.  Smith  and  Brown  [15]  applied  LES  with  a  Smagorinsky-type 
subgrid  closure  to  compute  the  complex  3-D  flow  patterns  around  building  arrangements 
in  a  2  km  x  2  km  section  of  Lower  Manhattan.  DeCroix  [16]  applied  LES  to  simulate  the 
wind  field  in  a  small  region  in  downtown  Salt  Lake  City. 

The  present  paper  considers  the  predictive  performance  of  two  high-Reynolds  number  k- 
c  models  when  applied  to  the  complex  developing  flow  over  and  within  an  array  of  3- 
D  buildings.  In  consequence,  comparisons  of  model  predictions  with  experimental  data 
obtained  from  a  very  detailed  wind  tunnel  experiment  will  be  undertaken,  with  primary 
emphasis  on  evaluation  of  the  predictive  accuracy  of  the  model  for  the  mean  velocity  and 
turbulence  upstream,  within,  and  downstream  of  the  3-D  building  array.  The  work  reported 
here  is  an  extension  of  an  earlier  effort  in  which  the  performance  of  variants  of  the  k-e  model 
with  both  linear  and  nonlinear  eddy-viscosity  formulations  was  investigated  for  the  flow 
within  and  over  an  array  of  2-D  buildings  [14].  Finally,  the  high-resolution  CFD  simulations 
undertaken  here  will  be  used  to  diagnose  the  dispersive  stresses  in  the  urban-type  roughness 
array.  The  importance  of  the  dispersive  stresses,  which  arise  from  the  spatial  inhomogeneity 
in  the  mean  flow,  in  relation  to  the  spatially-averaged  Reynolds  stresses  within  an  urban- 
type  roughness  array  will  be  assessed. 

High-resolution  RANS  modelling 


The  equations 

The  governing  equations  for  a  steady,  incompressible,  and  neutrally-stratified  (adiabatic) 
fluid  flow  based  on  the  Reynolds-averaged  Navier-Stokes  (RANS)  approach  are 


du, 

dxi 


0, 


(1) 


dujUi 

dxj 


dp  d2Uj  d  ,-j—f  \ 


dx 
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where  the  Reynolds  averaging  of  a  quantity  is  denoted  by  drawing  a  bar  over  the  quantity. 
Here,  ut  and  u'  are  the  mean  and  fluctuating  velocities  in  the  x^-direction,  respectively, 
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with  i  =  1,  2,  or  3  representing  the  streamwise  x,  spanwise  y,  or  vertical  z  directions; 
Xi  =  (x,  y,  z)  —  x:  Uj  —  (u,v,w);  Ui  =  ( u,v,w )  =  (Lr,  V.  TV);  ^  is  the  kinematic  viscosity 
and  p  is  the  kinematic  pressure. 

Reynolds-averaging  the  Navier-Stokes  equation  gives  rise  to  the  so-called  kinematic  Rey¬ 
nolds  stresses  which  are  defined  as  the  tensor  The  Reynolds  stresses  depend  on  the 

velocity  fluctuations  u(,  and  introduce  new  unknown  quantities  in  the  RANS  equations, 
and  therefore  these  equations  no  longer  constitute  a  closed  system.  In  order  to  close  the 
system  of  equations,  we  need  further  equations  describing  the  relationship  between  u[u'- 
and  the  quantities  Uj  and  p  that  we  seek  to  determine.  This  is  known  as  the  closure 
problem  in  turbulence  modelling.  One  of  the  simplest  turbulence  models  for  involves 
approximating  the  Reynolds  stress  components  by  analog}’  with  a  Newtonian  type  of  linear 
constitutive  relationship  between  the  turbulence  stress  and  the  mean  strain-rate  tensor. 
This  model  uses  the  Boussinesq  eddy- viscosity  approximation  given  by 


where  Vt  is  the  kinematic  eddy  viscosity,  k  =  \u\u\  is  the  turbulence  kinetic  energy  (TKE), 
and  Sij  is  the  Kronecker  delta  function. 


Dimensional  analysis  dictates  that  the  eddy  viscosity  i>t  be  determined  by  the  product  of 
a  turbulence  velocity  scale  and  a  turbulence  length  scale.  In  the  k-e  model  (where  e  is  the 
rate  at  which  TKE  is  converted  into  thermal  internal  energy),  the  turbulence  velocity  and 
length  scales  are  chosen  to  be  proportional  to  k1/2  and  k3^2/e,  respectively.  Hence,  in  the 
k-e  modelling  framework  at  high-Reynolds  (high-Re)  numbers,  the  eddy  viscosity  is  given 

by  2 

«  =  c£,  (4) 

where  is  a  closure  constant. 


The  modelled  transport  equations  for  the  turbulence  kinetic  energy  and  the  viscous  dissi¬ 
pation  rate  are  given  by 
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respectively.  Here.  P  is  the  production  of  k  defined  as 


(5) 

(6) 
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,  duj 
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Two  different  high-Re  number  k-e  models  are  considered  here,  which  differ  primarily  in 
the  manner  in  which  the  production  term  P  [viz.,  Equation  (7)]  is  approximated.  In  the 
standard  (STD)  k-e  model,  the  production  term  P  is  simply  modelled  as  [17] 


P  =  vt 
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where 


(9) 


S|  =  *(2%Sy)1/2,  Sij 


1  f  dui  duj\ 

2  v  dx  j  dxi ) 


Here,  5y  is  the  mean  rate-of-strain  tensor  (or,  equivalently,  the  symmetric  part  of  the 
mean  velocity  gradient  tensor).  The  Kato-Launder  (K-L)  k-e  model  [18]  was  designed 
to  circumvent  the  stagnation  point  anomaly.  The  standard  k-e  model  predicts  excessive 
levels  for  k  and  vt  near  the  stagnation  point  of  flows  impinging  on  walls  and  other  solid 
surfaces  (e.g.,  leading  edge  of  an  airfoil).  To  suppress  the  excessive  turbulence  energy 
production  in  the  stagnation  zone,  the  strategy  adopted  in  the  K-L  k-e  model  is  to  modify 
the  approximation  for  the  production  P  in  Equation  (8)  as  follows: 


P  =  CMe|S||fi|, 


(10) 


where 


e 


(2Qy  Qy)1/2,  Q  ij 


1  f  ditj  duj  \ 

2  \dxj  dxi )  ’ 


(11) 


with  O ^  being  the  mean  rate-of-rotation  (vorticity)  tensor  (or,  equivalently,  the  antisym¬ 
metric  part  of  the  mean  velocity  gradient  tensor)*  This  simple  modification  in  essence 
equates  TKE  production  to  zero  in  irrotational  flows.  Consequently,  in  zones  of  irrotational 
straining  such  as  near  a  stagnation  point  |0|  ~  0  which  eliminates  the  excessive  production 
of  k  here.  However,  it  should  be  borne  in  mind  that  the  Kato-Launder  modification  of 
P  introduces  an  inconsistency  between  the  approximation  of  the  Reynolds  stresses  in  the 
fc-equation  and  the  Reynolds-averaged  Navier-Stokes  equation  [the  latter  of  which  uses  the 
Boussinesq  linear  stress-strain  relationship  of  Equation  (3)  to  approximate  the  Reynolds 
stresses] . 


Together,  the  transport  equations  for  k  and  e  contain  five  closure  constants  (C^,  crei 
Cei,  and  C€2)  which  must  be  determined  before  the  equations  can  be  solved.  The  k-e  model 
(standard  and  Kato-Launder)  employs  values  for  the  closure  coefficients  that  have  been 
arrived  at  by  comprehensive  data  fitting  over  a  wide  range  of  canonical  turbulent  flows. 
The  closure  coefficients  for  the  k-e  model  are 

Cp  —  0.09,  ak  =  h  ac  =  1.3,  Cel  =  1.44,  Ce2  -  L92.  (12) 


A  closed- form  solution  for  the  k-e  model  can  be  obtained  for  the  neutral  wall  shear  layer. 
The  solution  gives 


U  —T~  log  2  +  ^j 

K>v 


(13) 


where  «*  =  (— u'u’')1^2  is  the  friction  velocity  and  B  is  a  constant  of  integration.  For  this 
analytical  solution,  we  find  an  implied  value  for  the  von  I<  arm  an  constant,  kV)  of 

k2v  =  s/c;{Ca-C*)ot.  (14) 

Using  the  closure  coefficient  values  for  the  k-e  model,  kv  assumes  a  value  of  0.43.  The 
experimental  values  for  kv  are  mostly  in  the  range  0.41  ±  0.2,  so  the  implied  value  of  kv  in 
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the  k-e  model  is  consistent  with  these  measurements.  The  stress-intensity  ratio,  ul/k,  is 
predicted  to  be  0.3  (using  =  0.09)  which  agrees  well  with  experimental  measurements 
of  this  quantity  in  many  shear  flows  [19]. 


Numerical  framework 


The  coupled,  nonlinear  system  of  six  partial  differential  equations  which  models  the  steady- 
state,  three-dimensional  turbulent  flow  field  was  solved  using  the  STREAM  code  [20].  The 
algorithm  implemented  in  this  code  employs  a  fully-collocated,  cell-centred  storage  arrange¬ 
ment  for  all  transported  properties.  Within  an  arbitrary,  but  structured,  non-orthogonal 
finite-volume  system,  the  velocity  vector  is  decomposed  into  its  Cartesian  components  to 
which  the  momentum  equations  relate.  Diffusive  volume-face  fluxes  are  discretized  using  a 
second-order  central  differencing  scheme.  The  higher-order  quadratic  upwind  interpolation 
for  convective  kinematics  (QUICK)  scheme,  described  by  Leonard  [21],  is  used  to  approx¬ 
imate  the  convective  volume-face  fluxes.  This  second-order  scheme  for  the  discretization 
of  the  net  convective  flux  through  a  control  volume  combines  the  second-order  accuracy  of 
a  central  differencing  scheme  with  the  stability  inherent  in  an  upwind  differencing  scheme 
by  using  in  each  direction  separately  a  quadratic  upstream  interpolation.  The  second-order 
accuracy  of  the  QUICK  scheme  minimizes  numerical  diffusion  errors  that  are  characteristic 
of  first-order  accurate  discretization  schemes  such  as  the  hybrid  differencing  scheme  or  the 
power-law  differencing  scheme  [22]. 


Mass  conservation  is  enforced  indirectly  by  solving  a  pressure-correction  equation  which,  as 
part  of  the  iterative  sequence,  steers  the  pressure  towards  a  state  in  which  all  mass  residuals 
in  the  cells  are  negligibly  small.  The  iterative  scheme  used  here  to  enforce  mass  conservation 
is  the  Semi-Implicit  Method  for  Pressure-Linked  Equations  (SIMPLE)  described  in  detail 
by  Patankar  [22].  The  SIMPLE  algorithm  provides  a  method  for  connecting  the  discretized 
forms  of  the  momentum  and  continuity  equations  to  give  an  equation  linking  the  pressure 
correction  at  a  node  to  its  neighboring  nodes.  In  this  iterative  solution  sequence,  Ui  ( i  = 
1, 2, 3)  are  initially  obtained  with  an  estimated  pressure  field.  This  is  continuously  updated 
by  reference  to  the  local  mass  residuals,  which  are  used  to  steer  the  pressure  field  towards 
the  correct  level.  Within  this  scheme,  the  transport  equations  for  U(  ( i  =  1,2,3),  k:  and 
e  and  the  pressure-correction  equation  are  solved  sequentially  and  iterated  to  convergence, 
defined  by  reference  to  LI  residual  norms  for  the  mass  and  momentum  components.  Here, 
the  LI  residual  norm  is  defined  as  the  sum  of  absolute  residuals  over  all  grid  points  of 
the  flow  domain.  The  residual  norm  provides  a  quantitative  measure  of  how  perfectly  the 
discretization  equations  are  satisfied  by  the  current  values  of  the  dependent  variables.  The 
LI  residual  norms  for  the  mass  and  momentum  components  were  normalized  by  the  mass 
and  momentum  fluxes  at  the  inflow  plane.  A  convergent  solution  was  assumed  after  each 
normalized  LI  residual  norm  decreased  below  0.001. 

A  fully  collocated  variable  storage,  in  combination  with  a  central  differencing  for  pressure, 
is  known  to  provoke  checkerboard  oscillations,  reflecting  the  pressure-velocity  decoupling. 
To  avoid  this,  the  widely  used  interpolation  practice  of  Rhie  and  Chow  [23]  is  adopted 
to  interpolate  cell-face  velocities  from  the  nodal  values.  This  method  for  interpolation 
essentially  introduces  a  fourth-order  smoothing  term  based  on  the  pressure,  and  prevents 
the  occurrence  of  spurious  pressure  modes. 
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The  wind  tunnel  experiment 


The  wind  tunnel  experiment  is  fully  described  in  Brown  et  al.  [24],  and  only  the  important 
details  of  the  experiment  will  be  presented  here.  The  experiments  were  conducted  in  an 
open-return  type  wind  tunnel  at  the  US  Environmental  Protection  Agency's  (EPA)  Fluid 
Modelling  Facility  (FMC).  This  wind  tunnel  has  a  working  test  section  of  length  18.3  m, 
width  3.7  m,  and  an  adjustable  roof  of  height  approximately  2.1  m  to  eliminate  streamwise 
pressure  gradients  and  allow  for  a  non-accelerating  free  stream  flow. 

The  wind  tunnel  experiment  considered  in  this  paper  simulated  a  neutrally-stratified  atmo¬ 
spheric  boundary  layer  flow  over  an  array  of  3-D  buildings.  A  schematic  of  the  3-D  building 
array  used  for  the  wind  tunnel  experiment  is  shown  in  Figure  1  (top).  This  array  consisted 
of  sharp-edged  cubes  with  a  characteristic  dimension  of  L  =  W  =  H  =  0.15  m,  where  L, 
W ,  and  H  are  the  length,  width,  and  height  of  the  obstacle.  A  total  of  77  cubes  was  placed 
in  an  aligned  array  consisting  of  7  rows  of  11  cubes.  The  array  filled  the  entire  spanwise 
dimension  (width)  of  the  wind  tunnel.  The  streamwise  and  spanwise  spacing  between  cubes 
was  if,  giving  a  frontal  area  index  of  Xp  =  0.25,  where  Xf  =  WH/(SX  +  L)(Sy  -I-  W)  (Sx 
and  Sy  denote  the  streamwise  face-to-face  spacing  and  spanwise  face-to-face  spacing  of  the 
obstacles,  respectively).  The  building  array  was  immersed  in  a  simulated  neutral  atmo¬ 
spheric  boundary  layer  which  was  created  in  the  tunnel  using  spires  and  floor  roughness 
with  a  roughness  length  of  approximately  0.001  m.  This  combination  produced  a  simulated 
boundary  layer  with  a  depth  of  1.8  m  and  a  friction  velocity  iu  ~  0.24  ms-1.  The  mean 
streamwise  velocity  in  the  upstream  approach  flow  can  be  approximated  by  the  following 
power-law  form 

f  =  (!f6'  (15> 

where  tin  —  3  ms-1  is  the  reference  velocity  of  the  upstream  flow  at  z  =  H.  The  Reynolds 
number,  based  on  the  building  height  H  and  the  reference  velocity  tin,  is  30,000.  This 
reference  Reynolds  number  is  larger  than  the  lower  limit  of  4,000  required  for  Reynolds 
number  independence  in  a  shear  flow  around  surface-mounted  obstacles  [25]. 

The  three  velocity  components  were  measured  with  a  pulsed-wire  anemometer  (PWA)  using 
a  pulsing  rate  of  10  Hz  and  an  averaging  time  of  120  seconds  at  each  measurement  loca¬ 
tion.  The  velocity  measurements  consisted  of  81  vertical  profiles  of  the  mean  velocity  and 
turbulence  velocity  variance  in  the  three  coordinate  directions,  and  the  turbulence  kinetic 
energy.  These  profiles  were  measured  in  the  vertical  centre  (symmetric)  plane  of  the  3-D 
building  array  (i.e.,  in  the  vertical  x-z  plane  at  y/H  =  0).  The  profiles  extended  from  3.5 H 
upstream  of  the  windward  face  of  the  first  row  of  cubes  to  7 .5 H  downstream  of  the  leeward 
face  of  the  last  row  of  cubes.  All  streamwise  distances  are  referenced  relative  to  the  location 
of  the  upstream  edge  of  the  first  row  of  cubes  at  x/H  =  0. 

Flow  domain  discretization  and  boundary  conditions 


For  the  simulations  to  be  reported,  the  model  domain  used  to  simulate  the  flow  within  and 
over  the  3-D  building  array  spanned  —5  <  x/H  <28  with  the  windward  face  of  the  first 
row  of  buildings  placed  at  x/H  =  0.  The  height  of  the  domain  is  8 H.  Simulations  were 
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conducted  also  with  a  domain  height  of  10#  and  12#,  but  we  found  that  these  alternative 
domain  heights  had  a  negligible  effect  on  the  solution.  Our  simulations  are  focused  on  flow 
changes  within  and  above  the  urban  canopy  (to  a  height  of  about  3#).  A  domain  depth  of 
8#  was  sufficiently  deep  to  provide  an  appropriate  domain  within  which  the  flow  changes 
near  the  surface  can  occur  without  being  moulded  by  the  boundary  conditions  imposed  at 
the  outer  boundary-layer. 

The  discretization  of  the  flow  domain  is  exhibited  in  Figure  1  (middle  and  bottom).  The 
geometry  possesses  a  two-fold  symmetry  in  the  spanwise  direction.  More  specifically,  since 
the  geometry  is  symmetric  about  the  vertical  x-z  centre  plane  at  y/H  =  0  and  the  vertical 
x-z  plane  at  y/H  =  1  (i.e.,  vertical  x-z  plane  centred  in  the  streamwise  oriented  street 
canyon),  only  one-half  of  the  full  geometry  in  the  spanwise  direction  need  to  be  included 
in  the  computational  domain,  so  0  <  y/H  <  1.  The  simulation  was  carried  out  in  a 
three-dimensional  Cartesian  framework. 

To  investigate  the  sensitivity  of  the  flow  solution  to  the  discretization  of  the  flow  domain, 
two  different  grids  were  used.  In  particular,  a  coarse  grid  of  135  x  22  x  40  cells  (in  the 
streamwise,  spanwise,  and  vertical  directions)  and  a  fine  grid  of  202  x  32  x  60  cells  (in  the 
streamwise,  spanwise,  and  vertical  directions)  were  employed  for  the  simulations.  The  grid 
arrangement  adopted  for  the  fine  mesh  is  shown  in  Figure  1  (middle  and  bottom),  which 
shows  the  calculation  grid  in  the  x-z  (middle)  and  x-y  (bottom)  planes.  The  grid  lines  were 
preferentially  concentrated  near  the  solid  surfaces  (ground,  rooftop,  and  building  walls) 
where  the  gradients  in  the  flow  properties  are  expected  to  be  greatest,  and  the  spacing 
between  the  grid  lines  was  gently  stretched  with  increasing  distance  from  the  solid  surfaces. 
In  particular,  the  grid  spacing  in  the  upstream  and  downstream  blocks  (viz.,  blocks  upstream 
of  the  windward  face  of  the  first  row  of  buildings  and  downstream  of  the  leeward  face  of  the 
last  row  of  buildings)  were  expanded  in  the  streamwise  direction  away  from  the  building 
faces  by  constant  factors  of  1.1  and  1.05,  respectively.  The  grid  spacing  for  the  block  above 
the  top  of  the  array  (viz.,  for  z/H  >  1)  was  stretched  in  the  vertical  direction  away  from 
the  top  of  the  buildings  by  a  constant  factor  of  1.1. 

The  dependence  of  the  flow  solutions  on  grid  resolution  was  investigated  using  the  coarse 
and  fine  grids.  Figure  2  shows  an  example  of  the  prediction  of  the  vertical  profiles  of  mean 
streamwise  velocity  U  and  turbulence  kinetic  energy  k  obtained  with  the  coarse  and  fine 
mesh.  These  profiles  were  calculated  at  x  =  1425  mm  ( x/H  =  9.5)  in  the  vertical  centre 
plane  (y/H  =  0)  of  the  array.  This  x-location  is  centred  in  the  spanwise  oriented  street 
canyon  between  the  5th  and  6th  rows  of  cubes.  The  agreement  between  the  coarse-  and 
fine-grid  solutions  is  very  good,  and  generally  results  for  U  were  found  to  display  negligible 
grid  sensitivity  for  the  two  sets  of  grids  identified  above.  However,  the  results  for  k  showed  a 
sensitivity  (albeit  relatively  small)  to  the  coarse-  and  fine- grid  resolutions.  Generally,  these 
results  imply  that  the  simulations  that  will  be  presented  in  the  following  section  using  the 
fine  mesh  can  be  regarded  as  essentially  grid-independent. 

Since  the  computational  domain  is  chosen  to  be  large  compared  with  the  3-D  array  of 
obstacles,  the  flow  at  its  boundary  is  far  removed  from  the  region  affected  by  the  obstacle 
array.  Hence,  free  boundary  conditions  are  imposed  at  all  air-to-air  boundaries  in  the  flow 
domain.  At  the  inflow  (inlet)  boundary,  the  measured  profiles  of  undisturbed  streamwise 
mean  velocity,  u  =  uj,  and  turbulence  kinetic  energy,  fc,  are  used.  The  vertical  mean 
velocity  w  =  =  0  at  the  inflow  plane  at  x/H  —  —5.  Unfortunately,  the  dissipation  rate  e 
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was  not  measured  and,  consequently,  was  estimated  by  assuming  an  equilibrium  turbulent 
flow  in  which  the  rates  of  production  and  destruction  of  turbulence  are  in  near  balance,  so 
e  «  C^4k3/2/(kvz)  where  kv  m  0.4  is  von  Karman’s  constant  and  z  is  the  normal  distance 
from  the  wall. 

At  the  lateral  (spanwise)  boundaries,  a  symmetry  boundary  condition  (i.e.,  no  flow  across 
the  boundary  and  no  scalar  flux  across  the  boundary)  was  imposed  along  the  vertical  x-z 
(symmetry)  planes  at  y/H  =  0  and  y/H  =  1  (cf.  Figure  1).  Far  downstream,  at  the  outflow 
boundary,  the  flow  was  assumed  to  reach  a  fully-developed  state  where  no  changes  occur 
in  the  flow  direction.  Hence,  at  the  outflow  boundary,  the  gradients  of  all  flow  variables 
are  assumed  to  be  zero  in  the  flow  direction;  viz.,  du/dx  =  dvfdx  =  dwj'dx  =  dk/dx  — 
defdx  =  0.  At  the  upper  boundary,  we  used  free-slip  conditions  for  all  flow  variables. 

At  all  the  solid  boundaries  (ground,  obstacle  walls,  obstacle  roofs),  standard  wall  functions 
are  applied  for  the  mean  velocities  and  turbulence  quantities.  Wall  functions  are  used 
to  reduce  the  computational  cost  associated  with  the  alternative  of  using  a  low-Reynolds 
number  turbulence  model  to  numerically  integrate  the  solution  through  the  entire  near¬ 
wall  region  including  the  viscous  sublayer  up  to  the  wall  where  no-slip  and  impermeability 
conditions  can  be  applied.  Typically  this  will  require  a  very  fine  grid  near  the  wall  because 
spatial  variations  in  near-wall  turbulence  structure  are  large  here  due  to  the  combined 
influence  of  viscosity  and  wall-induced  anisotropy  Log-law  based  ‘wall  laws'  are  adopted  in 
conjunction  with  high-Reynolds  number  turbulence  models  to  bridge  the  viscous  sublayer. 
Basically  the  flow  quantities  at  the  first  grid  point  above  the  solid  surface,  which  is  located 
outside  the  viscous  sublayer,  are  related  to  the  wall  friction  velocity  based  on  the  assumption 
of  a  semi-logarithmic  velocity  distribution,  A  local  equilibrium  of  turbulence  (where  the 
production  and  dissipation  of  turbulence  are  in  balance)  is  assumed  to  prevail  near  the  wall. 
When  these  flaw  of  the  wall’  type  boundary  conditions  are  applied,  the  diffusive  flux  of  k 
through  the  wTall  is  zero,  yielding  a  boundary  condition  that  the  normal  derivative  of  k  is 
zero  at  the  wall.  The  particular  specification  and  implementation  of  optimum  near-wall 
relationships,  for  use  with  the  k-e  model  (obtained  from  extensive  computing  trials),  have 
been  described  by  Versteeg  and  Malalasekera  [26] . 

Results  and  discussion 


Plots  of  the  mean  flow  streamlines  (curves  whose  tangent  is  parallel  to  the  mean  velocity 
vector  projected  into  the  given  plane  at  each  point  in  this  plane)  in  the  impingement  zone 
upstream  of  the  building  array  and  through  the  first  three  building  rows  are  displayed  in 
Figure  3  for  the  STD  and  K-L  k-e  models.  Similarly,  the  mean  streamline  patterns  of 
airflow  in  the  exit  region  downstream  of  the  last  row  of  buildings  obtained  with  the  STD 
and  K-L  k-e  models  are  exhibited  in  Figure  4,  These  mean  flow  patterns  were  obtained 
along  the  vertical  centre  plane  of  the  array  at  y/H  —  0.  In  addition,  isopleths  of  turbulence 
kinetic  energy  k  have  been  superimposed  on  the  mean  flow  streamlines  in  Figures  3  and 
4.  The  mean  flow  streamline  plots  are  particularly  useful  for  determining  the  shape  of 
the  recirculation  regions,  and  helpful  for  identifying  the  location  of  primary  and  secondary 
eddies  in  the  complex,  3-D  flow  considered  here. 

It  is  informative  to  contrast  these  results  for  the  3-D  building  array  with  those  for  the  2-D 
building  array  reported  in  our  previous  study  [14],  We  recall  that  the  2-D  building  array 


8 


DRDC  Suffield  TR  2004-108 


consisted  of  rectangular  blocks  of  equal  height  and  length  (H  =  L  =  0.15  m)  with  each 
block  extending  from  wall-to-wall  in  the  spanwise  direction  across  the  wind  tunnel.  In  the 
3-D  building  array,  the  high  levels  turbulence  energy  are  confined  in  the  region  extending 
from  the  windward  rooftop  edge  of  the  first  row  of  buildings  to  about  the  windward  rooftop 
edge  of  the  second  row  of  buildings  (see  Figure  3).  This  is  markedly  different  than  the 
2-D  building  array,  where  high  levels  of  turbulence  energy  were  observed  near  and  above 
the  top  of  the  array  resulting  in  a  ‘plume’  of  increased  turbulence  energy  extending  from 
the  windward  edge  of  the  first  building  to  about  the  leeward  edge  of  the  third  building 
in  the  array,  with  the  centreline  of  the  TKE  ‘plume’  at  z/H  «  1.25.  The  reason  for  this 
difference  arises  from  the  fact  that  the  displacement  of  streamlines  around  the  individual 
buildings  in  the  first  row  of  the  3-D  array  reduces  the  production  of  k  in  the  stagnation 
region  upstream  of  the  windward  face  of  the  array.  In  contrast,  in  the  2-D  building  array, 
the  lateral  displacement,  of  streamlines  is  not  possible  and  the  flow  in  the  impact  region  is 
necessarily  forced  over  the  first  building.  In  addition,  lateral  displacement  of  streamlines 
by  the  individual  buildings  in  the  3-D  array  also  suppresses  the  formation  of  a  separation 
zone  on  the  rooftops  of  the  buildings  in  the  first  row.  However,  the  deflection  of  streamlines 
upward  over  the  first  building  in  the  2-D  array  is  much  more  severe  than  in  the  3-D  array 
with  a  greater  streamline  compression  over  the  top  of  the  first  building.  Hence,  in  the  2-D 
array,  a  leading-edge  separation  bubble  on  the  roof  emananting  from  the  windward  rooftop 
edge  ( x/H  =  0)  of  the  first  building  is  evident,  with  a  reattachment  point  at  x/H  ~  0.8. 
The  prediction  of  these  qualitative  features  of  the  flow  over  and  through  the  2-D  and  3-D 
building  arrays  by  the  k-e  models  is  consistent  with  measurements  of  the  velocity  field  in 
these  two  arrays  conducted  by  Brown  et  al.  ([24],  [27]). 

The  standing  vortex  recirculation  pattern  in  the  street  canyons  of  the  2-D  array  is  markedly 
different  than  that  in  the  3-D  array.  In  the  2-D  array,  the  stagnation  point  (vortex  core) 
of  the  standing  vortex  is  approximately  located  at  half  canopy  height  (z/H  ss  0.5)  and 
half-way  between  two  consecutive  rows  of  buildings  in  the  spanwise  oriented  street  canyon. 
In  contrast,  for  the  3-D  array,  the  centre  of  the  recirculation  zone  between  two  consecutive 
rows  of  buildings  is  located  closer  towards  the  leeward  wall  of  the  upstream  building  and 
upwards  towards  the  rooftop  level  to  about  z/H  zz  0.75.  This  prediction  is  consistent  with 
velocity  vector  plots  of  the  vortex  recirculation  flow  in  the  3-D  array  presented  bv  Brown 
et  al.  [24], 

The  displacement  of  the  vortex  core  in  the  recirculation  that  develops  within  the  street 
canyons  of  the  3-D  array  appears  to  be  the  result  of  two  zones  of  swirling  flow  originating 
from  the  ground  surface  m  the  large  separation  region  that  develops  behind  a  cube  of  the 
array  as  exhibited  in  Figure  5.  Here,  a  pair  of  counter-rotating  vortices  in  the  horizontal 
x-y  plane  at  z/H  —  0.75  can  be  seen,  corresponding  to  two  areas  of  swirling  flow  which  are 
the  ‘footprints’  of  an  arch-shaped  vortex  that  forms  in  the  recirculating  flow  behind  a  3-D 
building.  This  arch  vortex  has  been  elucidated  in  flow  visualization  studies  and  detailed 
Laser  Doppler  Anemometry  measurements  of  the  velocity  field  around  an  isolated  cube  [28]. 
A  pictoral  depiction  of  the  complex  nature  of  the  flow  around  an  isolated  cube  is  given  in 
Hosker  [6]. 

The  topology  of  the  flow  pattern  of  the  recirculation  zone  behind  the  last  row  of  buildings 
in  the  3-D  array  is  different  than  that  in  the  2-D  array.  More  specifically,  the  trailing 
recirculation  zone  in  the  2-D  building  array  does  not  include  the  presence  of  a  spiral  saddle 
point  observed  in  the  trailing  recirculation  zone  of  the  3-D  building  array  (see,  Figure  4). 
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The  reader  is  referred  to  Hunt  et  al.  [29],  Tobak  and  Peake  [30],  and  Globus  et  al,  [31]  for 
further  information  on  the  application  of  topologically  based  kinematieal  concepts  to  the 
study  of  complex  flow  patterns.  In  addition,  the  recirculating  bubble  that  forms  downstream 
of  the  last  building  in  the  2-D  array  is  larger  than  that  in  the  3-D  array.  The  measured 
reattachment  point  occurs  at  x/H  ~  3.8  downstream  of  the  leeward  face  of  the  last  building 
in  the  2-D  array,  whereas  in  the  3-D  array  the  measured  reattachment  point  occurs  at 
x / H  ^  1,5  downstream  of  the  leeward  face  of  the  last  row  of  buildings.  Interestingly,  this 
measured  reattachment  length  for  the  recirculation  zone  behind  a  building  in  the  last  row  of 
the  3-D  array  is  almost  exactly  equal  to  that  measured  for  flow  over  an  isolated  cube  where 
the  reattachment  length  of  the  separation  zone  behind  the  cube  was  found  to  be  x/H  1.6 
[28].  Both  the  STD  and  K-L  k-e  models  predict  a  reattachment  length  (cf.  Figure  4)  for 
the  recirculation  zone  in  the  exit  region  of  x/H  ~  2.5,  returning  a  reattachment  length 
that  overpredicts  the  experimental  value  by  about  1.0 H.  Furthermore,  the  wind  shear  that 
develops  in  the  thin  layer  that  borders  the  separation  line  emanating  from  the  leeward 
rooftop  edge  of  the  last  building  in  the  2-D  array  is  larger  than  that  in  the  3-D  array.  As 
a  result  of  higher  shear  straining  here  and,  hence,  greater  production  of  turbulence  kinetic 
energy,  the  peak  values  of  k  in  the  exit  region  from  15  <  x/H  <  18  at  z/H  «  0.75  are 
larger  in  the  2-D  array  [14]  than  in  the  3-D  array. 

Finally.  Figures  3  and  4  provide  a  comparison  of  the  predicted  levels  of  the  turbulence 
energy  in  the  building  array  for  the  STD  and  K-L  k~e  models.  As  expected,  the  large  le\el 
of  production  of  turbulence  energy  near  the  stagnation  point  just  upstream  of  the  windward 
face  of  the  array  (impact  region)  predicted  by  the  STD  model  has  been  suppressed  in  the 
K-L  model.  In  the  exit  region,  however,  the  differences  in  the  levels  of  turbulence  energy 
predicted  by  the  STD  and  K-L  k-e  models  are  negligible. 

In  the  following  sub-sections,  the  vertical  profiles  of  the  horizontal  wind  component  L  (z), 
the  vertical  wind  component  IK(z),  and  the  turbulence  kinetic  energy  k  computed  within 
and  over  the  3-D  building  array  using  the  standard  and  Kato-Launder  fc-e  models  will  be 
compared  with  experimental  profiles  at  a  number  of  streamwise  locations.  These  longitudi¬ 
nal  locations  extend  from  the  impingement  zone  upstream  of  the  array,  through  the  array 
interior,  to  the  exit  region  downstream  of  the  array. 


Mean  streamwise  velocity 

The  vertical  profiles  of  mean  streamwise  velocity,  u  =  17,  in  the  vertical  centre  plane  at 
12  selected  ^-locations  upstream,  through,  and  downstream  of  the  3-D  building  array  are 
displayed  in  Figures  6-8.  Agreement  of  the  predictions  provided  by  the  STD  and  K-L 
jfc-e  models  with  the  experimental  data  is  excellent  at  almost  all  the  r-locations.  Both 
turbulence  closure  models  provide  virtually  identical  predictions  for  u. 

At  x  =  —75  mm  (x/H  =  —0.5)  in  the  impact  region  upstream  of  the  windward  array  face, 
the  wind  speed  reduction  below  and  the  wind  speed  increase  above  the  level  z/H  «  4/3 
relative  to  the  undisturbed  upstream  mean  wind  is  reproduced  well  (cf.  Figure  6).  The 
speed-up  of  the  flow  over  the  roof  of  a  building  in  the  first  row  [i.e,  at  x  =  75  mm  ( x/H 
0.5)]  and  the  absence  of  a  separation  zone  with  flow  reversal  on  the  roof  here  agrees  well 
with  the  experimental  data.  At  x  =  225  mm  ( x/H  —  1.5),  which  is  at  the  centre  of  the 
first  street  canyon,  the  magnitude  of  the  very  strong  wind  shear  at  or  near  canopy  height 
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and  the  magnitude  of  the  reverse  flow  in  the  canyon  below  about  z/H  «  2/3  are  correctly 
predicted  by  the  model  simulations. 

The  vertical  profiles  of  both  predicted  and  measured  u  shown  in  Figure  7  (for  the  array 
interior)  correspond  to  the  fully-developed  flow  within  and  over  the  array.  Note  that  the 
measured  (predicted)  u  at  x  =  1125  mm  (x/H  =  7.5)  in  the  fourth  street  canyon  is  almost 
the  same  as  that  measured  (predicted)  at  x  =  1425  mm  (x/H  =  9.5)  in  the  fifth  street 
canyon.  Also,  u  measured  and  predicted  at  x  =  1275  mm  ( x/H  =  8.5)  over  the  rooftop  of  a 
building  in  the  fifth  row  is  almost  the  same  as  that  measured  and  predicted  at  x  =  1575  mm 
(x/H  —  10.5)  over  the  rooftop  of  a  building  in  the  sixth  row.  In  particular,  the  measured 
u  in  the  fourth  street  canyon  and  over  the  rooftop  of  a  building  in  the  fifth  row  differ  by 
less  than  about  1%  (which  is  within  the  measurement  uncertainties)  from  the  measured 
u  in  a  comparable  position  in  the  fifth  street  canyon  and  over  the  rooftop  of  a  building 
in  the  sixth  row.  Consequently,  vertical  profiles  of  u  appear  to  be  fully-developed  (viz., 
longitudinal  modulations  in  u  from  one  street  canyon  to  the  next  appear  to  have  reached 
an  approximate  streamwise  equilibrium)  by  about  the  fourth  street  canyon,  with  the  result 
that  the  flow  in  the  array  interior  downstream  of  the  fourth  street  canyon  could  have  been 
modelled  by  simply  imposing  periodic  boundary  conditions  in  the  streamwise  direction  for  a 
flow  domain  that  consists  of  a  single  row  of  buildings  and  the  associated  downstream  street 
canyon. 

In  the  exit  region  downstream  of  the  leeward  face  of  the  array,  the  mean  streamwise  velocity 
below  the  level  z/H  &  2/3  is  underpredicted  by  the  flow  models  [cf.  Figure  8  where  it  is 
seen  that  the  models  predict  a  larger  reverse  mean  velocity  than  what  is  actually  observed 
at  x  =  2175  mm  (or,  x/H  =  14.5)  implying  an  elongation  of  the  recirculation  zone,  and 
underpredict  u  at  x  =  2475  mm  (or,  x/H  =  16.5)].  This  defect  may  be  the  result  of 
the  inability  of  the  models  to  capture  the  subtle  interactions  between  the  curvature  strain 
and  turbulent  stresses  in  the  curved  shear  layer  bordering  the  recirculation  zone  (i.e.,  an 
underprediction  of  the  shear  stress  in  this  layer  would  elongate  the  recirculation  zone).  The 
latter  discrepancy  provokes  another  weakness  in  the  solutions — namely,  an  insufficient  rate 
of  momentum  recovery  in  the  far  wake  region  following  reattachment. 


Mean  vertical  velocity 


The  vertical  profiles  of  mean  vertical  velocity,  w  =  W,  in  the  vertical  centre  plane  (y/H  =  0) 
at  12  selected  x-locations  upstream,  through,  and  downstream  of  the  3-D  building  array 
are  displayed  in  Figures  9-11.  The  predicted  results  generally  agree  fairly  well  with  the 
measurements  in  a  qualitative  sense,  but  the  conformance  between  predictions  and  mea¬ 
surements  is  not  as  good  as  those  for  the  mean  streamwise  velocity. 

As  the  front  face  of  a  building  in  the  first  row  is  approached,  the  vertical  velocity  is  positive 
as  air  is  forced  up  and  over  the  rooftop  of  the  building.  This  qualitative  feature  in  the 
observations  is  reproduced  quite  well  by  the  modelled  results,  as  is  the  vertical  position  of 
the  maximum  positive  mean  vertical  velocity  at  or  near  the  top  of  the  building  (z/H  1)  [cf. 
measured  and  predicted  mean  vertical  velocities  at  x  —  —75  mm  (x/H  —  —0.5)  in  Figure  9]. 
However,  the  magnitude  of  the  maximum  positive  vertical  velocity  is  overpredicted  by  about 
25%  here.  Indeed,  upstream  of  the  windward  face  of  a  building  in  the  first  row,  the  model 
predicts  that  there  is  a  vertical  pressure  gradient  with  a  high  pressure  region  centred  at 
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about  z/H  «  2/3  that  forces  the  air  over  the  rooftop  of  the  building  above  this  level 
and  down  towards  the  ground  below  this  level.  We  note  that  the  streamline  curvature 
over  a  building  rooftop  in  the  first  row  as  well  as  immediately  upstream  of  this  rooftop  is 
concave  upwards  (cf.  Figure  3),  implying  a  destabilizing  (exchange-enhancing)  influence  on 
the  turbulence. 


The  upward  motion  of  air  above  the  rooftop  of  a  building  in  the  first  row  predicted  by  the 
models  is  in  conformance  with  the  wind  tunnel  results,  but  the  magnitude  of  the  positive 
mean  vertical  velocity  is  overpredicted  here  by  about  20%  [cf.  the  vertical  profiles  of  W 
measured  and  predicted  at  x  =  75  mm  (x/H  =  0.5)].  The  negative  vertical  velocity  close 
to  the  rooftop  surface  is  captured  in  the  numerical  simulation  (albeit,  negative  W  here  is 
underpredicted).  The  prediction  of  the  downward  (negative)  vertical  velocity  at  x  =  225 
mm  (or,  x/H  =  1.5  which  is  at  the  centre  of  the  first  street  canyon)  is  in  good  conformance 
with  the  observations. 


Vertical  profiles  of  computed  and  measured  mean  vertical  velocity  W  within  the  fourth  street 
canyon  and  over  the  rooftop  of  a  building  in  the  fifth  row  are  identical  (approximately  or 
better)  to  the  profiles  obtained,  respectively,  within  the  fifth  street  canyon  and  rooftop  of 
a  building  in  the  sixth  row.  Hence,  consistent  with  the  behaviour  of  [/,  the  mean  vertical 
velocity  appears  to  be  fully  developed  by  the  fourth  street  canyon  in  the  array.  Figure  10 
exhibits  W  in  the  array  interior  where  the  mean  flow  is  fully  developed. 


Figure  11  compares  computed  and  measured  profiles  of  IT  in  the  last  street  canyon  and  in 
the  exit  region  of  the  array.  Generally  speaking,  computed  and  measured  W  in  the  exit 
region  (x/H  >  13.0)  are  largely  reproduced  by  the  model  predictions  with  the  discrepancies 
between  the  predictions  and  observations  being  less  than  about  5%,  except  within  the 
recirculation  bubble  in  the  near  wake  region.  For  example,  at  x  =  21  (5  mm  ( x/H  =  14.5,  or 
1.5#  downstream  of  the  leeward  face  of  the  array),  the  magnitude  of  the  downward  vertical 
velocity  at  z/H  «  1  is  underpredicted  by  about  60%.  However,  the  prediction  of  the  weak 
upward  motion  below  about  z/H  ~  1/3  here  is  in  good  conformance  with  the  measurements. 
The  rate  of  recovery  of  the  mean  vertical  momentum  beyond  the  reattachment  point  in  the 
far  wake  region  to  its  far  upstream  reference  state  is  reproduced  accurately. 


Turbulence  kinetic  energy 


Vertical  profiles  of  turbulence  kinetic  energy,  A,  at  the  same  x-locations  as  those  for  the 
mean  streamwise  and  vertical  velocities  are  displayed  in  Figures  12  to  14.  Solutions  for 
turbulence  energy  derived  from  the  standard  and  Kato-Launder  fe-e  models  differ  primarily 
in  the  impact  region  upstream  of  the  windward  face  of  the  array  and  over  the  rooftop  of 
a  building  in  the  first  row.  In  particular,  consistent  with  the  results  exhibited  earlier  in 
Figure  3,  the  fc-profiles  in  Figure  12  at  x  =  —75  mm  ( x/H  —  —0.5)  and  at  x  —  75  mm 
(x/H  —  0.5)  show  lower  turbulence  energy  levels  for  the  K-L  model  predictions  than  for 
the  STD  model  predictions.  The  excessive  turbulence  kinetic  energy  just  upstream  of  a 
building  in  the  first  row  is  absent  in  the  K-L  model  results,  so  there  is  less  turbulence 
energy  swept  downstream  into  the  region  over  the  rooftop  of  this  building.  In  consequence, 
at  both  ^-stations,  the  K-L  model  predictions  are  seen  to  be  in  better  conformance  with 
the  measurements  than  the  STD  model  predictions. 
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Although  the  predicted  u-  and  ^-profiles  in  the  street  canyons  and  over  the  building  rooftops 
of  the  array  agree  quite  well  with  the  experimental  data,  the  conformance  between  the  pre¬ 
dicted  A:-profiles  and  the  measurements  is  not  so  satisfactory.  In  general  a  weakness  shared 
by  both  the  STD  and  K-L  k-e  models  is  that  they  underpredict  the  levels  of  turbulence 
energy  in  the  street  canyons  and  at  the  rooftops  of  the  buildings  in  the  array  interior.  The 
peak  value  of  k  [which  occurs  at  or  near  the  canopy  top  (z/H  &  1)]  decreases  monotonically 
in  the  streamwise  direction  from  the  first  street  canyon  and  reaches  a  near  constant  value 
by  the  about  the  fourth  street  canyon.  Here,  the  predicted  peak  value  of  k  is  about  65%  of 
the  measured  value. 

This  is  in  conformance  with  the  expected  performance  of  the  standard  k-e  model  for  pre¬ 
diction  of  a  boundary-layer  flow  over  a  flat  plate  as  reported  by  Patel  et  al.  [32]  in  their 
review  article.  Here,  it  was  remarked  that  the  peak  value  of  k  predicted  by  the  standard  k-e 
model  for  a  fully-developed  flow  over  a  flat  plate  is  about  70%  of  the  experimentally  mea¬ 
sured  value.  Patel  et  al.  [32]  suggest  that  this  deficiency  is  associated  with  the  fundamental 
weakness  of  the  linear  (Boussinesq)  stress-strain  relationship  related  to  the  single-scale 
assumption  and  the  omission  of  any  anisotropic  eddy-viscosity  effects.  In  principle,  this 
deficiency  can  be  circumvented  by  either  using  a  Reynolds  stress  transport  model  (Launder 
et  al.,  [33])  or,  as  an  effective  compromise  between  the  full  second-moment  closure  and  the 
two-equation  model  using  an  algebraic  stress  model  [34]  or  Durbin's  k-e-v2  model  [35].  In 
the  latter  model,  vt  in  Equation  (4)  is  replaced  by 

vt  =  C^v2,  (16) 

with  v2  (defined  to  be  the  Reynolds  stress  component  normal  to  the  mean  flow  streamline) 
obtained  from  its  own  transport  equation.  The  essential  purpose  of  Durbin’s  k-e-v2  model 
is  to  represent  the  tendency  of  a  solid  surface  to  suppress  the  turbulent  fluctuations  in 
the  direction  normal  to  the  surface  without  having  to  resort  to  the  full  complexity  of  the 
Reynolds  stress  transport  equations. 

The  streamw'ise  evolution  of  k  is  dominated  by  the  vertical  spreading  of  the  strong  shear 
layer  that  develops  at  or  near  the  top  of  the  urban  canopy  (z/H  «  1).  At  the  centre  of  this 
layer  lying  just  above  the  building  height  H,  the  peak  values  of  k  attenuate  downstream, 
with  the  turbulence  energy  at  the  central  core  of  the  layer  being  “exported"  to  regions  below 
and  above  the  core  by  turbulent  and/or  pressure  diffusion,  causing  the  turbulence  energy 
to  increase  in  these  regions. 

The  shape  of  the  A’q:>ro files  and  the  position  of  maximum  k  in  the  exit  region  [i.e.,  at  x  =  2175 
and  2475  mm  (or,  x/H  =  14.5  and  16.5)  in  Figure  14]  appear  to  be  correctly  predicted. 
The  quantitative  agreement  of  turbulence  energy  levels  with  the  experimental  data  is  also 
fairly  good  here  [e.g.,  the  peak  value  of  k  at  x  =  2475  mm  (x/H  =  16.5)  is  underestimated 
only  by  about  25%].  Also,  the  vertical  fc-profile  at  x  =  2475  mm  (x/H  =  16.5)  has  the 
appearance  of  a  detached  mixing  layer  in  spite  of  the  fact  that  the  mean  flow  has  reattached 
at  x/H  ~  1.5  downstream  of  the  leeward  face  of  the  array.  This  profile  form  is  a  residue 
of  the  strong  shear  layer  that  detaches  from  a  building  in  the  last  row  of  the  array  and 
spreads  vertically  outward  by  pressure  and  turbulent  diffusion  until  this  “mixing”  layer 
makes  contact  with  the  ground  surface. 
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Diagnosis  of  dispersive  stresses 


After  validation,  the  model  can  be  used  to  diagnose  and  to  evaluate  the  importance  of  the 
dispersive  stresses  relative  to  the  usual  spatially-averaged  Reynolds  stresses  within  and  over 
the  urban-type  roughness  array  studied  here.  The  dispersive  stresses  arise  from  the  spatial 
inhomogeneity  of  the  mean  velocity  field  in  the  building  array*  and  can  be  inferred  directly 
from  the  predicted  flow  quantities  obtained  from  the  model  simulations. 

The  available  experimental  evidence  suggests  that  in  fine-scaled  plant  canopies*  the  dis¬ 
persive  stress  is  negligible  in  comparison  to  the  spatially-averaged  Reynolds  stresses  near 
and  above  the  canopy  top  (Raupach  et  al.,  [36]).  However*  Bohm  et  al.  [37]  recently 
demonstrated  that  in  a  model  plant  canopy,  the  dispersive  stress  is  comparable  to  the 
spatially-averaged  Reynolds  stresses  near  the  bottom  of  the  canopy  (z/H  <  0.5),  but  both 
stresses  were  found  to  be  small  here  owing  to  the  decay  of  these  stresses  with  depth  from  the 
top  of  the  canopy.  Cheng  and  Castro  [38]  measured  the  dispersive  stress  above  urban-like 
roughness  arrays  (coarse-scaled  cuboid  arrays)  and  found  them  to  be  negligible  compared 
to  the  spatially-averaged  Reynolds  stresses  for  z/H  >  1,  which  is  consistent  with  the  re¬ 
sults  of  Raupach  et  al  [36]  and  Bohm  et  al  [37]  for  plant  canopies.  However,  owing  to 
experimental  difficulties,  Cheng  and  Castro  were  not  able  to  measure  the  dispersive  stresses 
within  an  urban-type  canopy  and  concluded  that  whether  the  dispersive  stress  is  very  small 
compared  to  the  spatially- averaged  Reynolds  stress  here  “remains  an  open  question'5  that 
“requires  very  extensive  and  technically  difficult  measurements  before  a  conclusive  answer 
is  possible55 . 


Let  us  define  the  spatial  average  of  a  flow  quantity  <j>  as 


v 


x  +  r,  i)  dr, 


where  V  =  Vf  +  Vs  is  the  averaging  volume  consisting  of  both  a  fluid  subvolume  Vf  and  a 
solid  subvolume  14*  In  the  following,  V  will  be  taken  to  be  a  thin  horizontal  slab  (chosen  to 
resolve  the  characteristic  variation  of  the  flow  property  <j>  in  the  vertical)  that  encompasses 
one  or  more  obstacles  in  the  array.  Given  the  spatial  averaging  operation  in  Equation  (17), 
the  flow  quantity  <f>  can  be  decomposed  into  its  spatial  average  and  a  departure  therefrom 
to  give 

<t>  =  (<£)+ f,  (18) 

where  the  departure  satisfies  {o")  —  0.  In  &  multiply  connected  space  (e.g.,  urban  canopy), 
spatial  differentiation  and  spatial  averaging  do  not  commute,  so 


dtp 

dxi 


d(4>) 

dxi 


(19) 


where  S  is  the  part  of  the  averaging  volume  V  that  coincides  with  the  obstacle  surfaces  and 
rii  is  the  unit  normal  vector  pointing  from  V  into  S  (Raupach  et  al.,  [36]). 


If  we  apply  the  spatial  averaging  operator  of  Equation  (17)  to  the  Reynolds-averaged  Navier- 
Stokes  equation  [cf.  Equation  (2)]  and  use  Equation  (18),  we  get 


d(Uj)(Ui) 

dxj 


d(p) 


(20) 
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where 


(21) 


and 


viscous  drag  form  drag 


(22) 


The  form  (pressure)  and  viscous  drag  terms,  defined  in  Equation  (22),  arise  directly  as  a 
consequence  of  the  non-commutation  of  spatial  differentiation  and  spatial  averaging.  The 
spatially-averaged  total  stress  tensor  in  Equation  (21)  includes  a  spatially-averaged  kine¬ 
matic  Reynolds  stress  (u-u')  and  a  kinematic  dispersive  stress  (u”u")  that  arises  from  the 
spatial  correlations  in  the  mean  velocity  field  caused  by  horizontal  heterogeneity  of  the 
canopy. 


In  this  section,  the  kinematic  dispersive  normal  and  shear  stresses,  (uf/uf/)  and  will 

be  diagnosed  for  a  coarse-scaled  array  of  cubes  (cf.  Figure  1)  using  the  high-resolution  RANS 
results  presented  in  Section  5  and  their  importance  relative  to  the  associated  kinematic 
spatially-averaged  Reynolds  stresses  (u'u;)  and  ( ufw ')  will  be  evaluated.  To  this  purpose, 
the  3-D  building  array  described  earlier  is  first  decomposed  into  seven  ‘averaging  units’,  each 
unit  having  a  streamwise  length  Lx  =  2 H  that  encompasses  a  row  of  buildings  in  the  array 
and  the  associated  downstream  spanwise  oriented  street  canyon  as  shown  in  Figure  15(a). 
Each  ‘averaging  unit’  is  further  decomposed  into  a  number  thin  horizontal  slabs  of  constant 
volume  V  that  are  stacked  in  the  vertical  (or  z-)  direction  as  shown  in  Figure  15(b)  (each  of 
which  constitutes  an  averaging  volume  V  that  will  be  used  to  define  the  spatial  averages  for 
the  flow  quantities  in  the  array).  Hence,  by  reference  to  Figure  15(b),  the  spatial  average 
of  a  (time-averaged)  flow  variable  <j>  will  be  calculated  explicitly  as 

j  ry~H  nx—xo^2H 

(<I>)(Z)  =  7?7H  /  </>dxdy  (23) 

J  y— 0  J  x=Xq 

with  the  understanding  that  <j)  =  0  within  the  building.  In  Equation  (23),  periodicity  of 
the  flow  in  the  array  in  the  spanwise  direction  has  been  used.  Finally,  xo  =  2 (i  —  1  )H 
(i  =  1,2, ...  ,7)  for  ‘averaging  unit’  #i. 


First,  let  us  consider  the  mean  streamwise  and  vertical  velocities,  <j>  =  u  and  w,  respec¬ 
tively.  There  are  720  vertical  profiles  of  u  and  w  calculated  in  each  ‘averaging  unit’  for  the 
simulations  reported  earlier.  As  an  example,  the  vertical  profiles  of  u  and  iu  and  their  corre¬ 
sponding  spatial  averages  (u)  and  (w)  for  averaging  unit  #5  are  exhibited  in  Figures  16(a) 
and  (c),  respectively.  The  flow  in  the  array  is  fully  developed  in  averaging  unit  #5  (ar¬ 
ray  interior),  implying  that  d(u)/dx  =  0  and  d(w)/dz  —  0  (by  the  continuity  equation), 
so  ( w)(z )  =  0  (since  (w)(z  =  0)  =  0  by  the  impermeability  condition  at  the  ground  sur¬ 
face).  The  expected  vanishing  of  (w)  in  averaging  unit  #5  (where  the  canopy  flow  is  fully 
developed)  is  consistent  with  the  result  shown  in  Figure  16(c). 

Figures  16(b)  and  (d)  show  vertical  profiles  of  the  fluctuating  velocities  (i.e.,  departures  of 
the  time-averaged  velocities  from  their  spatial  averages)  u"  and  w"  and  their  corresponding 
spatial  averages  (u")  and  (w"),  respectively.  As  seen,  both  (u")  and  (w")  are  zero,  which 
is  consistent  with  the  fact  that  for  the  choice  of  the  averaging  volume  V  in  Equation  (15), 
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(ptlsj  =  o.  From  the  results  displayed  in  Figures  16(b)  and  (d),  vertical  profiles  of  0  u  and 
unwn  can  be  constructed  and  are  given  in  Figures  16(e)  and  (f),  respectively*  which  include 
also  the  spatial  averages  of  these  profiles  over  the  Averaging  unit3*  The  latter  quantities 
correspond  to  the  kinematic  dispersive  normal  stress,  {ut,uH),  and  shear  stress.  { uNw,f ), 
respectively. 


Comparison  of  vertical  profiles  of  the  kinematic  dispersive  normal  stress,  (u'V').  and  the 
spatial  average  of  the  kinematic  Reynolds  normal  stress,  (uV),  for  Averaging  units3  #1  to 
#6  in  the  building  array  are  shown  in  Figures  17-19.  Consistent  with  Cheng  and  Castro  s 
experimental  observations  [38],  {uffuff)  is  negligibly  small  compared  to  {ufuf)  above  the 
canopy  height  (z/H  >  1).  Below  the  canopy  height,  however,  (u,fuff)  is  approximately 
3-6  times  larger  than  implying  that  the  dispersive  stresses  are  significant  within  an 

urban-type  roughness  array  (e.g.,  the  coarse-scaled  cuboid  array  considered  here).  However, 
this  comparison  assumes  that  the  predictions  of  vfu*  provided  by  the  k-e  model  (which  is 
based  on  the  Boussinesq  eddy- viscosity  hypothesis)  are  correct.  The  k-e  model  was  shown 
earlier  to  under  predict  the  turbulence  energy  levels,  by  as  much  as  a  factor  of  two.  implying 
that  ufu*  could  be  underestimated  by  as  much  as  a  factor  of  two.  Even  so,  the  results  of 
Figures  17-19  would  still  imply  that  (unuN)  is  comparable  in  magnitude  to  ( ufu ')  within 
this  urban-type  roughness  array. 


Vertical  profiles  of  the  spatially-averaged  Reynolds  shear  stress  (ufwf)  and  the  dispersive 
shear  stress  { uf,w,f )  for  the  first  6  Averaging  units3  in  the  array  are  displayed  in  Figure  20- 
22.  The  results  here  clearly  demonstrate  that  although  the  dispersive  shear  stress  is  very 
small  compared  to  the  spatially-averaged  Reynolds  shear  stress  at  and  above  the  canopy 
top  (consistent  with  measurements  made  by  Cheng  and  Castro  [38]  in  this  region),  the 
magnitude  of  the  dispersive  shear  stress  is  comparable  to  that  of  the  spatially-averaged 
Reynolds  shear  stress  below  the  canopy  top.  Again,  it  should  be  emphasized  that  this 
conclusion  implicitly  assumes  that  the  values  of  ufwf  predicted  by  the  k-e  model  are  correct 
for  this  building  array.  Unfortunately,  no  data  for  the  shear  stress  are  available  from  the 
wind  tunnel  experiment  to  test  this  assumption.  However,  the  use  of  the  linear  (Boussinesq) 
relationship  between  stresses  and  strains  is  adequate  for  the  prediction  of  shear  stresses  in 
thin  shear  flows,  and  even  if  \v/n7\  is  underpredicted  by  a  factor  of  two  in  our  case,  the 
conclusion  here  would  still  be  valid. 


Interestingly,  we  found  that  (uNw,t)  and  (ufwf)  are  of  opposite  sign  for  z/H  <  1.  and  the 
magnitude  of  {uHwft)  is  larger  than  that  of  ( utwt )  for  the  first  three  Averaging  units  (i.e., 
units  #1  to  #3).  For  the  succeeding  three  Averaging  units3  (i.e.,  units  #4  to  #6),  the 
magnitude  of  both  stresses  are  comparable.  Furthermore,  the  shape  of  |{u/tt,/)j  is  different 
than  that  of  \{uf,wf%  This  is  an  interesting  and  important  observation,  because  it  is 
not  (' v!9w” ),  that  directly  affects  the  momentum  flux.  As  an  example,  let  us  focus 
on  (u,wf)  and  ( u,fwff )  diagnosed  for  Averaging  unit3  ^5  where  the  flow  is  fully-developed 
(approximately  or  better).  Note  that  here  the  vertical  gradients  of  these  two  stresses  for 
0  <  z/H  <  0.5  are  of  opposite  sign  and  of  similar  magnitude,  implying  that  §z(uf,w,t)  is 
'destructive3  relative  to  (u¥)  (viz.,  the  momentum  flux  arising  from  the  turbulent  shear 
stress  is  reduced  by  the  dispersive  shear  stress).  On  the  other  hand,  for  0.5  <  z/H  <  1, 
-g^{uuwtf)  has  the  same  sign  as  u*wr ),  and  the  dispersive  shear  stress  here  is  ‘constructive 

in  nature  (i.e.,  the  momentum  flux  arising  from  the  turbulent  and  dispersive  shear  stresses 
reinforce  each  other). 
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Conclusions 


Two  variants  of  k-e  turbulence  closure  schemes  (namely,  the  standard  and  Kato-Launder 
models)  were  used  to  simulate  the  flow  over  a  3-D  building  array,  with  a  plan  and  frontal  area 
density  of  0.25,  corresponding  to  the  ‘skimming  flow*  regime.  The  predictions  from  these  two 
models  were  compared  with  available  data  obtained  from  a  wind  tunnel  experiment.  The 
agreement  between  the  predicted  mean  streamwise  velocity  profiles  and  the  experimental 
data  is  generally  very  good,  with  the  greatest  discrepancy  in  u  occurring  in  the  recirculation 
zone  immediately  downstream  of  the  leeward  face  of  the  array.  Qualitatative  features  in 
the  mean  vertical  velocity  have  been  reproduced  well,  and  the  quantitative  agreement  of 
predicted  w  with  the  experimental  data  is  also  fairly  good  throughout  most  of  the  flow 
domain,  with  the  largest  discrepancy  between  measured  and  predicted  values  of  w  being 
about  25%. 

The  peak  values  of  the  turbulence  kinetic  energy  within  the  building  array  (i.e.,  over  the 
rooftops  of  the  buildings  and  the  tops  of  the  street  canyons  in  the  array)  are  generally 
under  predicted.  Although  the  turbulence  kinetic  energy  does  not  match  the  experimental 
data  to  quite  the  same  degree  as  the  mean  streamwise  and  vertical  velocities,  we  note  that 
the  quantitative  agreement  between  predicted  and  measured  k  is  within  a  factor  of  two  at 
worst,  but  more  usually  within  about  50%.  Fortunately,  the  prediction  of  the  mean  velocity 
field  (17,  V.  W)  appears  to  be  relatively  insensitive  to  the  inaccuracies  in  the  predicted  value 
of  k. 

The  comparison  of  predicted  and  observed  flow  quantities  in  the  building  array  indicates 
that  a  simple  RANS  model  (k-e  turbulence  closure)  of  the  complex  mean  flow  through  a 
regular  building  array  is  accurate  enough  to  warrant  its  application  to  the  prediction  of 
mean  flows  within  and  over  other  more  complex  urban-type  roughness  arrays.  However, 
if  a  more  accurate  prediction  of  turbulence  quantities  is  needed,  then  more  sophisticated 
RANS  models  based  on  the  transport  of  Reynolds  stresses  may  be  required.  It  is  expected 
that  the  latter  class  of  models  would  be  able  to  predict  the  effects  of  streamline  curvature, 
separation,  recirculation,  and  reattachment  on  normal  stress  anisotropy  better,  albeit  at  a 
substantially  increased  level  of  mathematical  and  computational  complexity. 

The  high-resolution  RANS  simulation  results,  after  careful  validation  against  the  wind  tun¬ 
nel  experiment,  were  subsequently  used  to  diagnose  the  dispersive  stresses  within  and  over 
the  building  array.  In  conformance  with  Cheng  and  Castro’s  [38]  experimental  observations 
over  urban-type  roughness  arrays,  the  dispersive  stresses  were  found  to  be  negligibly  small 
(i.e.,  insignificant)  in  comparison  to  the  associated  spatially-averaged  Reynolds  stresses 
above  the  canopy  height  (i.e.,  z/H  >  1).  However,  although  Cheng  and  Castro  were  unable 
to  measure  the  dispersive  stresses  within  their  urban-type  roughness  arrays,  we  have  applied 
our  numerical  simulations  to  evaluate  the  importance  of  dispersive  stresses  within  a  building 
array  consisting  of  coarse-scaled  cuboids.  It  was  found  that  the  dispersive  stresses  (arising 
from  the  spatial  inhomogeneity  in  the  mean  flow)  are  comparable  to  the  spatially-averaged 
Reynolds  stresses  within  the  building  array  (i.e.,  z/H  <  1).  Interestingly,  the  vertical  gra¬ 
dients  of  the  dispersive  shear  stress  in  the  array  interior  where  the  flow  is  fully  developed 
have  opposite  and  similar  parity  to  the  vertical  gradients  of  the  spatially-averaged  Reynolds 
shear  stress,  respectively,  for  z/H  <  0.5  and  z/H  >  0.5.  This  implies  that  the  kinematic 
momentum  fluxes  contributed  by  the  dispersive  shear  stress  act  destructively  relative  to 
that  contributed  by  the  spatially-averaged  Reynolds  shear  stress  in  the  bottom  half  of  the 
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canopy,  whereas  in  the  top  half  both  shear  stresses  act  constructively  to  determine  the  total 
momentum  fluxes. 

Given  that  dispersive  stresses  are  comparable  to  spatially-averaged  Reynolds  stresses  within 
an  urban-type  roughness  array  (coarse-scaled  urban  canopy),  the  problem  that  needs  to  be 
confronted  is  the  construction  of  a  physically-based  model  for  the  dispersive  stress  tensor 
Tf.  =  (u"u")  that  can  be  used  to  predict  the  spatially-averaged  time-mean  velocity  in  the 
canopy  [^Equations  (20)  and  (21)].  A  number  of  important,  rigorous  properties  wjll  need 
to  be  taken  into  account  in  specifying  an  appropriate  model  for  rfj.  Certainly,  ry  needs 
to  be  positive  semi-definite  (viz.,  rfj  needs  to  be  realizable,  so  that  the  dispersive  kinetic 
energy  kd  =  | rd  is  a  positive  quantity).  Because  the  dynamics  represented  by  Equation  (20) 
is  Galilean  invariant,  the  dispersive  stress  gradients  drfj/dxj  which,  physically,  are  forces 
acting  on  the  spatially-averaged  time-mean  flow  need  to  be  Galilean  invariant  also  (and, 
hence,  a  model  for  these  dispersive  stress  gradients  should  at  least  maintain  the  Galilean 
invariance).  The  dispersive  stress  tensor  must  depend  explicitly  on  the  characteristic  filter 
width  A  =  Vl/Z  [cf.  Equation  (17)]  used  to  spatially  average  the  flow  (and,  certainly,  we 
expect  t? .  -•>  0  as  A  -»  0+).  Finally,  the  modelling  of  the  dispersive  stresses  is  complicated 
by  the  fact  that  rfj  can  locally  be  either  diffusive  (dissipative)  in  nature  [implying  an  energy 
drain  from  the  resolved  scales  to  the  sub- filter  scales  (i.e.,  scales  smaller  than  the  filter  width 
A)]  or  reactive  (dispersive)  in  nature  [implying  a  backward  scatter  of  energy  from  the  sub¬ 
filter  scales  to  the  resolved  scales].  In  particular,  note  for  the  specific  case  shown  in  Figure  22 
that  rf3  =  ( u"w ")  has  an  opposite  sign  to  ( u'w ')  so  that  rf3d{u)/ dz  >  0.^  In  other  words, 
the  shear  production  term  which  converts  resolved-scale  kinetic  energy  §  (ui){ui)  to  sub¬ 
filter  scale  kinetic  energy  ufuf)  would  correspond  here  to  a  backscatter  of  energy  from 
the  sub-filter  scales  to  the  resolved  scales  (i.e.,  scales  greater  than  A). 
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Figure  1.  A  three-dimensional  perspective  view  (top)  showing  the  geometry  of  the  aligned 

array  of  cubes  (7  rows  of  cubes  with  11  cubes  in  each  row)  and  the  Cartesian  coordinate  * 

system  used.  Here,  x  is  in  the  streamwise  direction,  y  is  in  the  spanwise  direction,  and  z  is 

in  the  vertical  direction.  The  x-z  view  (middle)  at  y  —  0  and  x-y  view  (bottom)  at  z  =  0 

of  the  computational  grid  used  for  the  prediction  of  flow  over  the  aligned  array  of  cubes  are  * 

shown.  This  Cartesian  grid  arrangement  corresponds  to  the  202  x  32  x  60  (fine)  grid. 
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Figure  2.  Grid-sensitivity  analysis:  vertical  profiles  of  the  mean  streamwise  velocity  U 
and  the  turbulence  kinetic  energy  k  obtained  at  x  =  1425  mm  (x/H  =  9.5)  between  5th 
and  6th  rows  of  cubes  in  the  array  for  the  coarse  and  fine  grids. 


DRDC  Suffield  TR  2004-108 


23 


k 

0.0384 

0.0355 

0.0326 

0.0297 

0.0269 

0.0240 

0.0211 

0.0182 

0.0154 

0.0125 

0.0096 

0.0067 

0.0038 

0.0010 


3 

2.5 
2 

1.5 
1 

0.5 

0 

-0.5 


k 

0.0384 

0.0355 

0.0326 

0.0297 

0.0269 

0.0240 

0.0211 

0.0182 

0.0154 

0.0125 

0.0096 

0.0067 

0.0038 

0.0010 


Figure  3.  Mean  flow  streamlines  in  the  impingement  zone  upstream  of  and  through  the 
first  three  rows  of  buildings  of  the  array  obtained  with  the  standard  and  Kato-Launder  k-e 
turbulence  closure  models.  Turbulence  kinetic  energy  isopleths  have  been  superimposed 
on  the  mean  streamline  pattern  of  airflow.  The  flow  patterns  and  k  were  obtained  in  the 
vertical  centre  plane  of  the  array  at  y/H  =  0. 


* 


24 


DRDC  Suffield  TR  2004-108 


Figure  4.  Mean  flow  streamlines  in  the  exit  region  downstream  of  the  last  row  of  buildings 
of  the  array  obtained  with  the  standard  and  Kato-Launder  k-e  turbulence  closure  models. 
Turbulence  kinetic  energy  isopleths  have  been  superimposed  on  the  mean  streamline  pattern 
of  airflow.  The  flow  patterns  and  k  were  obtained  in  the  vertical  centre  plane  of  the  array 
at  y/H  =  0. 
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Figure  5.  Velocity  vector  plot  in  the  horizontal  x-y  plane  at  z/H  —  0,75  obtained  in  the 
vicinity  of  a  building  in  the  6th  row  of  the  3-D  building  array.  This  plot  shows  two  areas 
of  swirling  flow  near  the  two  side  edges  of  the  building  which  are  the  footprints  of  an  arch 
vortex. 
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Figure  6.  Vertical  profiles  of  the  mean  streamwise  velocity  u  =  U  at  four  x-locations 
(x  =  — 75,  75,  225  and  375  mm),  obtained  from  a  high-resolution  numerical  simulation 
using  two  different  k-c  turbulence  closure  models,  are  compared  with  time-averaged  wind 
tunnel  measurements  at  the  same  locations. 
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Figure  7.  Vertical  profiles  of  the  mean  streamwise  velocity  u  =  U  at  four  allocations 
(x  =  1125,  1275,  1425  and  1575  mm),  obtained  from  a  high-resolution  numerical  simulation 
using  two  different  k-e  turbulence  closure  models,  are  compared  with  time-averaged  wind 
tunnel  measurements  at  the  same  locations. 
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Figure  8.  Vertical  profiles  of  the  mean  streamwise  velocity  u  =  U  at  four  ^-locations 
(x  =  1725,  1875,  2175  and  2475  mm),  obtained  from  a  high-resolution  numerical  simulation 
using  two  different  k-e  turbulence  closure  models,  are  compared  with  time-averaged  wind 
tunnel  measurements  at  the  same  locations. 
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Figure  9.  Vertical  profiles  of  the  mean  vertical  velocity  w  =  W  at  four  ^-locations  (a: 
—75,  75,  225  and  375  mm),  obtained  from  a  high-resolution  numerical  simulation  using 
two  different  k-t  turbulence  closure  models,  are  compared  with  time-averaged  wind  tunnel 
measurements  at  the  same  locations. 
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Figure  10.  Vertical  profiles  of  the  mean  vertical  velocity  w  =  W  at  four  ^-locations 
(x  =  1125,  1275,  1425  and  1575  mm),  obtained  from  a  high-resolution  numerical  simulation 
using  two  different  k-e  turbulence  closure  models,  are  compared  with  time-averaged  wind 
tunnel  measurements  at  the  same  locations. 
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Figure  11.  Vertical  profiles  of  the  mean  vertical  velocity  w  =  W  at  four  a>locations 
(x  =  1725,  1875,  2175  and  2475  mm),  obtained  from  a  high-resolution  numerical  simulation 
using  two  different  k-e  turbulence  closure  models,  are  compared  with  time-averaged  wind 
tunnel  measurements  at  the  same  locations. 
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Figure  12.  Vertical  profiles  of  the  turbulence  kinetic  energy  k  at  four  x-locations  (x  = 
—75,  75,  225  and  375  mm),  obtained  from  a  high-resolution  numerical  simulation  using 
two  different  k-e  turbulence  closure  models,  are  compared  with  time-averaged  wind  tunnel 
measurements  at  the  same  locations. 
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Figure  13.  Vertical  profiles  of  the  turbulence  kinetic  energy  k  at  four  x-locations  (x  =  1125 
1275,  1425  and  1575  mm),  obtained  from  a  high-resolution  numerical  simulation  usinf 
two  different  k-e  turbulence  closure  models,  are  compared  with  time-averaged  wind  tunne 
measurements  at  the  same  locations. 
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Figure  14.  Vertical  profiles  of  the  turbulence  kinetic  energy  k  at  four  x-locations  (x  =  1725, 
1875,  2175  and  2475  mm),  obtained  from  a  high-resolution  numerical  simulation  using 
two  different  k-e  turbulence  closure  models,  are  compared  with  time-averaged  wind  tunnel 
measurements  at  the  same  locations. 
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Figure  15.  Decomposition  of  the  building  array  into  (a)  seven  ‘averaging  units’,  each  unit 
of  which  consists  of  a  row  of  buildings  and  the  associated  downstream  spanwise  oriented 
street  canyon  (top);  (b)  an  ‘averaging  unit’  is  further  decomposed  into  thin  horizontal  slabs 
at  a  number  of  different  levels  z  (middle  and  bottom). 
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Figure  16.  Vertical  profiles  of  (a)  u,  ( u ),  and  (b)  u," ,  (u")  for  ‘averaging  unit’  #5  within 
and  over  the  building  array.  Other  parts  of  this  figure  are  continued  on  the  next  page. 
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Figure  16.  Vertical  profiles  of  (c)  w,  (w),  and  (d)  w",  (w")  for  ‘averaging  unit’  #5  within 
and  over  the  building  array.  Other  parts  of  this  figure  are  continued  on  the  next  page. 
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Figure  16.  Vertical  profiles  of  (e)  u"v!\  ( u" u //),  and  (f)  u"w " ,  (u"w")  for  ‘averaging  unit’ 
#5  within  and  over  the  building  array. 
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Figure  17.  Vertical  profiles  of  the  spatially-averaged  Reynolds  normal  stress  (u'u')  and 
the  dispersive  normal  stress  (u"u")  in  ‘averaging  units’  #1  and  #2  in  the  building  array. 
Further  parts  of  this  figure  for  other  ‘averaging  units’  are  continued  on  the  next  page. 
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Figure  18.  Vertical  profiles  of  the  spatially-averaged  Reynolds  normal  stress  (i u'u ')  and 
the  dispersive  normal  stress  {u"u")  in  ‘averaging  units’  #3  and  #4  in  the  building  array. 
Further  parts  of  this  figure  for  other  ‘averaging  units’  are  continued  on  the  next  page. 
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Figure  19.  Vertical  profiles  of  the  spatially-averaged  Reynolds  normal  stress  (u'u1)  an 
the  dispersive  normal  stress  (u"u")  in  ‘averaging  units’  #5  and  #6  in  the  building  array. 
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Figure  20.  Vertical  profiles  of  the  spatially-averaged  Reynolds  normal  stress  (mV)  and 
the  dispersive  normal  stress  ( u"u> ")  in  ‘averaging  units’  #1  to  #6  in  the  building  array. 
Further  parts  of  this  figure  for  other  ‘averaging  units’  are  continued  on  the  next  page. 
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Figure  21.  Vertical  profiles  of  the  spatially-averaged  Reynolds  normal  stress  (u'w1)  and 
the  dispersive  normal  stress  ( u"w ")  in  ‘averaging  units’  #3  to  #4  in  the  building  array. 
Further  parts  of  this  figure  for  other  ‘averaging  units’  are  continued  on  the  next  page. 


Ad 


DRDC  Suffield  TR  2004-108 


Figure  22.  Vertical  profiles  of  the  spatially-averaged  Reynolds  normal  stress  (u'w')  and 
the  dispersive  normal  stress  (u"w")  in  ‘averaging  units’  #5  to  #6  in  the  building  array. 
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